#include "utils.h"
/// Factory for symbol creation
const symbol & get_symbol(const string & s)
{
  static std::map<string, symbol> directory;
  std::map<string, symbol>::iterator i = directory.find(s);
  if (i != directory.end())
    return i->second;
  else
    return directory.insert(make_pair(s, symbol(s))).first->second;
}

lst set2lst(const exset & s) {
  lst l;
  exset::const_iterator it = s.begin();
  while (it != s.end()) {
    l.append(*it);
    ++it;
  }
  return l;
}

exset lst2set(const lst & l) {
  exset s;
  lst::const_iterator it = l.begin();
  while (it != l.end()) {
    s.insert(*it);
    ++it;
  }
  return s;
}

int loops_in_1pi_grpah(Graph& g) // 1PI graph => only one connected component(K) L=E-N+K
{
  return num_edges(g)-num_vertices(g)+1;
}

/// FIND PARALLEL EDGES (NOT USED)
bool has_par_edges(Graph::edge_descriptor edis,Graph G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  Graph::edge_iterator ei,ei_end;
  int par = 0;
  for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
    if((G[*ei].to_vertex == G[edis].to_vertex && G[*ei].from_vertex == G[edis].from_vertex) ||(G[*ei].to_vertex == G[edis].from_vertex && G[*ei].from_vertex == G[edis].to_vertex)) par ++;
  if (par > 1)return true;
  else return false;
}

// return vector of u-v parallel edges;
std::vector<Graph::edge_descriptor> edges_u_v(const Graph::vertex_descriptor& u,const Graph::vertex_descriptor& v,const Graph& G)
{
  std::vector<Graph::edge_descriptor> out_vec;
  Graph::out_edge_iterator pari,pari_end;
  Graph::edge_descriptor edis;
  std::vector<Graph::edge_descriptor> par_edge_set;
  for(tie(pari,pari_end)=out_edges(u,G);pari != pari_end;++pari)
    if(target(*pari,G) ==v)
      out_vec.push_back(*pari);
  return out_vec;
}

bool isNotConnected(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
	int n_v = num_vertices(G);
	std::vector<int> component(num_vertices(G));
	int num = connected_components(G, &component[0]);

	std::vector<int>::size_type i;
	//std::cout << "Total number of components: " << num << endl;
	//for (i = 0; i != component.size(); ++i)
          //	std::cout << "Vertex " << i <<" is in component " << component[i] << endl;
	std::cout << endl;
	if (num > 1) return true;
	else return false;
}
bool isWeakConnected(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  property_map < Graph, edge_component_t >::type
    wcomponent = get(edge_component, G);
  
  std::size_t num_comps = biconnected_components(G, wcomponent);
  

  if (num_comps > 1) 
    {
      // cout<<"Biconnected: "<<num_comps<<endl;
      return true;
    }
  else return false;
}

bool isConnected(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  
  if (isNotConnected(G) || isWeakConnected(G))return false;
  else return true;
}

int diagIndex(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  //write_graphviz(std::cout, G);
  //boost::property_map<Graph, edge_weight_t>::type PropogatorType = get(edge_weight, G);
  Graph Gp = G.root();
  typedef std::pair<Graph::edge_iterator, Graph::edge_iterator> edlst;
  Graph::edge_iterator f_edge,l_edge;
  Graph::vertex_iterator vi,vi_end;
  tie(f_edge,l_edge) = edges(G);
  int buf = 0;
  //Number of spinor lines calculation
  for (Graph::edge_iterator it = f_edge; it != l_edge; ++it)
    {
      //		buf += PropogatorType[*it];
      ///Only spinor lines summation t(G) = -1
      if(G[*it].type == 1) buf += G[*it].type;
      //		std::cout<<PropogatorType[*it]<<" ";      
      //      std::cout<<G[*it].type<<" ";
    }
  for (tie(vi,vi_end)=vertices(G); vi != vi_end; ++vi)
    {
      //		buf += PropogatorType[*it];
      ///Only spinor lines summation t(G) = -1
      if(G[*vi].type == 1) buf += G[*vi].type;
      //		std::cout<<PropogatorType[*it]<<" ";      
      std::cout<<"vtype"<<G[*vi].type<<" ";
    }
    int omega = 2*num_edges(G)-4*num_vertices(G)+buf+4;
    //int omega = -2*num_edges(G)+buf+4;
  //  std::cout<<std::endl<<"*********************************"<<std::endl;
  //std::cout<<"*** OMEGA = 2*L - 4*V + s + 4 = 2*"<<num_edges(G)<<" - 4*"<<num_vertices(G)<<" + "<<buf<<" + 4 = "<<omega<<std::endl;
  //std::cout<<"******************"<<omega<<"***************"<<std::endl;
  return omega;
}

/// if vertex is local or global external
bool is_external_vert(Graph::vertex_descriptor vdis,Graph& Gtmp)
{
  if(out_degree(vdis,Gtmp))
    {
      if(Gtmp.is_root())return(Gtmp[vdis].external);
      else  return((out_degree(vdis,Gtmp) < out_degree( Gtmp.local_to_global(vdis), Gtmp.root()))|| Gtmp.root()[Gtmp.local_to_global(vdis)].external);
    }
  else return false;
}

/// Check if vertex is REALLY end of spinor line
bool is_ext_spin_end(Graph::vertex_descriptor vdis,Graph& Gtmp)
{
  Graph::out_edge_iterator oei,oei_end;
  Graph::edge_descriptor edis;
  int n_sp = 0;
  for(tie(oei,oei_end) = out_edges(vertex(vdis,Gtmp),Gtmp);oei!=oei_end;++oei)
    if(Gtmp[*oei].type == 1)
      {
        n_sp++;
        edis = *oei;
      }
  cout<<" n_sp++ "<<n_sp<<endl;
  return  (1 == n_sp) && (vdis == vertex(Gtmp[edis].to_vertex,Gtmp) )&& is_external_vert(vdis,Gtmp) ;    
  // version with source and  target
  //  return ( (1 == n_sp) && (vdis == target(edis,Gtmp) ||vdis == source(edis,Gtmp)))&& is_external_vert(vdis,Gtmp) ;    
}

int num_ext_vert(Graph G)
{
  int n_ext = 0;
  Graph::vertex_iterator viu,viu_end;
  Graph::vertex_descriptor vdtmp;
  for (tie(viu, viu_end) = vertices(G); viu != viu_end; ++viu)
    if(is_external_vert(vdtmp=*viu,G) == true)n_ext++;
  return n_ext;
}


void print(Graph Gtmp)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  Graph::vertex_iterator viu,viu_end;
  Graph::vertex_descriptor vdtmp;
  cout<<"*********************************************"<<endl;
  cout<<"     PRINTING THE GRAPH ONLY FOR DEBUG!!!    "<<endl;
  cout<<"*********************************************"<<endl;
  print_graph(Gtmp);
  cout<< "Vert params: "<<endl;
  for (tie(viu, viu_end) = vertices(Gtmp); viu != viu_end; ++viu) 
    {
      cout<<setw(15)<<"External"<<is_external_vert(vdtmp = *viu,Gtmp)<<endl;
      cout<<setw(15)<<"External"<<Gtmp[*viu].visited<<endl;
      if(!Gtmp.is_root())cout<<setw(15)<<"Global dis"<<Gtmp.root()[Gtmp.local_to_global(*viu)].visited<<endl;
      if(!Gtmp.is_root())cout<<setw(15)<<"Momentum"<<Gtmp.root()[Gtmp.local_to_global(*viu)].gmomentum<<endl;
      else cout<<setw(15)<<"Momentum"<<Gtmp[*viu].gmomentum<<endl;
      cout<<setw(15)<<"Expression"<<Gtmp[*viu].gexpr<<endl;
      cout<<setw(15)<<"Local dis"<<*viu<<endl;
      if(!Gtmp.is_root())cout<<setw(15)<<"Global dis"<<Gtmp.local_to_global(*viu)<<endl;
      cout<<setw(17)<<"Num vert"<< num_vertices(Gtmp)<<endl;
      cout<<setw(17)<<"Num ext vert"<< num_ext_vert(Gtmp)<<endl;
      if(!Gtmp.is_root())cout<<setw(17)<<"Num vert in ROOT"<< num_vertices(Gtmp.root())<<endl;
 if(!Gtmp.is_root())cout<<setw(17)<<"Num ext vert in ROOT"<< num_ext_vert(Gtmp.root())<<endl;
      cout<<"---------------------------------------------"<<endl<<endl;
      //      cout<<setw(15)<<""<<Gtmp[*viu].<<endl;
    }
  cout<<"Edge params"<<endl;
  Graph::edge_iterator ei,ei_end;
  for(tie(ei, ei_end) = edges(Gtmp); ei != ei_end; ++ei)
    {
      cout<<setw(15)<<"From"<<Gtmp[*ei].from_vertex<<endl;
      cout<<setw(15)<<"To"<<Gtmp[*ei].to_vertex<<endl;
      cout<<setw(15)<<"Momentum"<<Gtmp[*ei].gmomentum<<endl;
      cout<<setw(15)<<"Expr"<<Gtmp[*ei].gexpr<<endl;
      cout<<setw(15)<<"Type"<<Gtmp[*ei].type<<endl;
      cout<<"---------------------------------------------"<<endl<<endl;
    }
  cout<<"*********************************************"<<endl;
}


ex ex_from_graph(Graph G)
{
 

  ///////////////////////////////////////////////
  ///  _____________________________________  ///
  /// |                                     | ///
  /// | Expression from graph, not subgraph | ///
  /// |_____________________________________| ///
  ///                                         ///
  ///////////////////////////////////////////////

  if(G.is_root())
    {
      GiNaC::numeric one1 = 1;
      ex out_expr = one1;
      cout<<"root gr Testing out_expr:"<< out_expr<<endl;
      // add out spinor lines and vertices with external legs 
      Graph::vertex_iterator vi,vi_end;
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) G[*vi].visited = false;
      Graph::edge_iterator ei,ei_end;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei)
        {
G[*ei].visited = false;
 cout<<"to "<<G[*ei].to_vertex<<" target"<<target(*ei,G)<<endl;
        }

      Graph::vertex_descriptor vdis;
      Graph::edge_descriptor edis;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && is_ext_spin_end(vertex(G[*ei].to_vertex,G),G) && !G[*ei].visited)
          {
            out_expr *= G[vertex(G[*ei].to_vertex,G)].gexpr;
            //    cout<<"01Vert in linesTesting out_expr:"<< out_expr<<endl;
            G[vertex(G[*ei].to_vertex,G)].visited = true;
            out_expr *= G[*ei].gexpr;
            // cout<<"01Edges  in lines Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
            edis = *ei;
            bool visited_spinor_line = false;
            while(!visited_spinor_line)
              {
                out_expr *= G[vertex(G[edis].from_vertex,G)].gexpr;
                //  cout<<"Vert in lines Testing out_expr:"<< out_expr<<endl;
                G[vertex(G[edis].from_vertex,G)].visited = true;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                bool lng = false;
                for(tie(out_ei,out_ei_end) = out_edges(vertex(G[edis].from_vertex,G),G);out_ei != out_ei_end;++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited)) 
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }

                if(lng == true)
                  {
                    out_expr *= G[*out_ei_tmp].gexpr;
                    //  cout<<"Edges in lines Testing out_expr:"<< out_expr<<endl;
                    G[*out_ei_tmp].visited = true;
                    edis = *out_ei_tmp;
                  }
                else visited_spinor_line = true;
              }
          }
      // spinor edges in loops
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && !G[*ei].visited)
          {
            bool visited_loop = false;
            edis = *ei;
            while(!visited_loop)
              {
                out_expr *= G[vertex(G[edis].to_vertex,G)].gexpr;
                //cout<<"Vertieces in loops Testing out_expr:"<< out_expr<<endl;
                G[vertex(G[edis].to_vertex,G)].visited = true;
                out_expr *= G[edis].gexpr;
                // cout<<"Edges inloops Testing out_expr:"<< out_expr<<endl;
                G[edis].visited = true;

                bool lng = false;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                for(tie(out_ei,out_ei_end) = out_edges(vertex(G[edis].from_vertex,G),G); out_ei != out_ei_end; ++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited))
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)edis = *out_ei_tmp;		
                else visited_loop = true;
	      
              }
          }
      // other edges and other verticies
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if(!G[*ei].visited)
          {
            out_expr *= G[*ei].gexpr;
            // cout<<"Other_edges Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
          }
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        if(!G[*vi].visited)
          {
            out_expr *= G[*vi].gexpr;
            // cout<<*vi<< "Other Vertices Testing out_expr:"<< out_expr<<endl;
            G[*vi].visited = true;
          }

      cout<<"Expr extracted from GRAPH no."<<get_property(G,graph_num)<<endl;
      //      cout<<"Not simplified: "<<out_expr<<endl;
      //cout<<"N = "<<out_expr.numer()<<endl;
      //cout<<"D = "<<out_expr.denom()<<endl;
      return out_expr;
    }
  ///////////////////////////////////////////////
  ///  _____________________________________  ///
  /// |                                     | ///
  /// | Expression from SUBgraph, not graph | ///
  /// |_____________________________________| ///
  ///                                         ///
  ///////////////////////////////////////////////
  else
    {
      //      print(G);
      GiNaC::numeric one1 = 1;
      ex out_expr = one1;
      cout<<"subgr Testing out_expr:"<< out_expr<<endl;
      // add out spinor lines and vertices with external legs 
      Graph::vertex_iterator vi,vi_end;
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) G[*vi].visited = false;
      Graph::edge_iterator ei,ei_end;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) G[*ei].visited = false;
        print(G);
        cout<<"After printing"<<endl;
      Graph::vertex_descriptor vdis;
      Graph::edge_descriptor edis;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && is_ext_spin_end(target(*ei,G),G) && !G[*ei].visited)
          {
            out_expr *= G[target(*ei,G)].gexpr;
             cout<<"01Vert in linesTesting out_expr:"<< out_expr<<endl;
             G[target(*ei,G)].visited = true;
            out_expr *= G[*ei].gexpr;
             cout<<"01Edges  in lines Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
            edis = *ei;
            bool visited_spinor_line = false;
            while(!visited_spinor_line)
              {
                out_expr *= G[source(edis,G)].gexpr;
                  cout<<"Vert in lines Testing out_expr:"<< out_expr<<endl;
                G[source(edis,G)].visited = true;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                bool lng = false;
                for(tie(out_ei,out_ei_end) = out_edges(source(edis,G),G);out_ei != out_ei_end;++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited)) 
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)
                  {
                    out_expr *= G[*out_ei_tmp].gexpr;
                      cout<<"Edges in lines Testing out_expr:"<< out_expr<<endl;
                    G[*out_ei_tmp].visited = true;
                    edis = *out_ei_tmp;
                  }
                else visited_spinor_line = true;
              }
          }
      // spinor edges in loops
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && !G[*ei].visited)
          {
            bool visited_loop = false;
            edis = *ei;
            while(!visited_loop)
              {
                out_expr *= G[target(edis,G)].gexpr;
                //   cout<<"Vertieces in loops Testing out_expr:"<< out_expr<<endl;
                G[target(edis,G)].visited = true;
                out_expr *= G[edis].gexpr;
                //   cout<<"Edges inloops Testing out_expr:"<< out_expr<<endl;
                G[edis].visited = true;

                bool lng = false;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                for(tie(out_ei,out_ei_end) = out_edges(source(edis,G),G); out_ei != out_ei_end; ++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited))
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)edis = *out_ei_tmp;		
                else visited_loop = true;
	      
              }
          }
      // other edges and other verticies
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if(!G[*ei].visited)
          {
            out_expr *= G[*ei].gexpr;
            //   cout<<"Other_edges Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
          }
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        if(!G[*vi].visited)
          {
            out_expr *= G[*vi].gexpr;
            cout<<*vi<< G[*vi].gexpr<< " Other Vertices Testing out_expr:"<< out_expr<<endl;
            G[*vi].visited = true;
          }

      cout<<"Expr extracted from GRAPH no."<<get_property(G,graph_num)<<endl;
      //      cout<<"Not simplified: "<<out_expr<<endl;
      //cout<<"N = "<<out_expr.numer()<<endl;
      //cout<<"D = "<<out_expr.denom()<<endl;
      return out_expr;

    }
}


// Some internal procedurs for permutation generation(e.g. increment in binary vector of booleans)
void inc(std::vector<bool>& vec)
{
  std::vector<bool>::iterator it;
  it = find(vec.begin(),vec.end(),false);
  *it=true;
  fill(vec.begin(),it,false);
}

bool is_full(std::vector<bool> vec)
{
  if (find(vec.begin(),vec.end(),false) == vec.end()) return true; 
  else return false;
}


Graphref_lst make_subgraphs(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |  Generating all divergent subgraphs | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  //num_vertices(G);
  cout<<endl<<"  MAKE_SUBGRAPHS PROCEDURE STARTED"<<endl<<endl;
  cout<<"gnum"<<get_property(G,graph_num)<<endl;
  cout<<ex_from_graph(G)<<endl;
  std::vector<bool> tmk(num_vertices(G),false);
  
  Graphref_lst l1;
  //tmk++;
  inc(tmk);
  
  while(!is_full(tmk))
    {
      //inc(tmk);  
      Graph::vertex_iterator vi,vi_end,viu,viu_end;
      Graph& Gtmp = G.create_subgraph();
      
      
      for (tie(viu, viu_end) = vertices(G); viu != viu_end; ++viu) 
	{
	  if(tmk[*viu])
	    {
	      add_vertex(vertex(*viu,G),Gtmp);
              Graph::vertex_descriptor vd = *viu;
              //              if(is_external_vert(vd,Gtmp))cout<<setw(18)<<"  "<<*viu<<"is ext"<<endl;
              //	      cout<<"Vertex "<<*viu<<" added to GRAPH with expr:";//<<Gtmp[*viu].gexpr<<endl;
	    }
	  //else cout<<"Vertex "<<*viu<<" not added"<<endl;
        }
      //M(Gtmp);
      
      
      if(isConnected(Gtmp)&&(diagIndex(Gtmp) >= 0)&&(num_edges(Gtmp) > 1))
        {
          /// Add External flag to vertices
          // cout<<std::right<< "Befor marking "<<setw(6)<<num_ext_vert(Gtmp)<<"from "<<num_vertices(Gtmp)<<" <<<<";
          // mark_ext_vert(Gtmp);
          // cout<<"After marking "<<setw(6)<<num_ext_vert(Gtmp)<<"from "<<num_vertices(Gtmp)<<" <<<<"<<endl;
          // print(Gtmp);
          //print_graph(Gtmp);
          //          cout<<setw(13)<<"Number of EDGES: "<<num_edges(Gtmp)<<endl;
          cout<<ex_from_graph(Gtmp)<<endl;
          get_property(Gtmp,graph_num) = get_property(G,graph_num);
          cout<<"original graph:"<<endl;
          print_vertices(G,get(vertex_index,G));
          //print_edges(G,get(edge_index,G));
          print_graph(G);
          cout<<"M(G):"<<endl;
          print_vertices(Gtmp,get(vertex_index,Gtmp));
          //print_edges(Gtmp,get(edge_index,Gtmp));
          print_graph(Gtmp);
          //cout<<" End of M, start T(g)"<<endl;
          //T(Gtmp);
          cout<<"sgnum"<<get_property(Gtmp,graph_num)<<endl;
          print_graph(Gtmp);
          //          Graph::edge_descriptor edis = *(edges(G).first);
          l1.push_back(ref(Gtmp));
        }

      inc(tmk);
    }

  /// IMPORTANT SWITCH BACK!!!!!
  l1.reverse();
  cout<<"Found "<<l1.size()<<" divergent subgraphs "<<endl;
  return l1;
  
  //Graph& G1 = G.create_subgraph();
}
 
/// END of 1PI subgraphs generation functions


bool isSubgraph(Graph& g1,Graph& g2)  // True if G2 is subgraph of G1
{
  std::vector<int> v1,v2;
	Graph::vertex_iterator vi,vi_end;
	for (tie(vi, vi_end) = vertices(g1); vi != vi_end; ++vi) 
	  v1.push_back(g1.local_to_global(*vi));
	std::sort(v1.begin(),v1.end());
	for (tie(vi, vi_end) = vertices(g2); vi != vi_end; ++vi) 
	  v2.push_back(g2.local_to_global(*vi));	
	std::sort(v2.begin(),v2.end());
	
	return std::includes(v1.begin(),v1.end(),v2.begin(),v2.end());
}

/// checks graph intersection
bool isIntersected(Graph& g1,Graph& g2)  // True if G2 is intersected with G1 or G1and G2 has common parts
{
  std::vector<int> v1,v2;
  std::vector<int>::iterator it_begin,it_end;
  Graph::vertex_iterator vi,vi_end;
  for (tie(vi, vi_end) = vertices(g1); vi != vi_end; ++vi) 
    v1.push_back(g1.local_to_global(*vi));
  std::sort(v1.begin(),v1.end());
  for (tie(vi, vi_end) = vertices(g2); vi != vi_end; ++vi) 
    v2.push_back(g2.local_to_global(*vi));	
  std::sort(v2.begin(),v2.end());
  std::vector<int> v3(v1.size()+v2.size());
  it_end = std::set_intersection(v1.begin(),v1.end(),v2.begin(),v2.end(),v3.begin());
  if (it_end == v3.begin()) return false;
  else return true;
  
}

bool isStrongIntersected(Graph& G1,Graph& G2)
{
  return (isIntersected(G1,G2) && !(isSubgraph(G1,G2) || isSubgraph(G2,G1)));
}


/// From list like G1,G2,G3...Gn
/// creates list of terms like G1*G2,G1*G3...
sorted_sum_of_M_prod make_terms(Graphref_lst l1)  
{
  cout<<num_vertices(unwrap_ref(l1.front()))<<endl;
  std::vector<bool> tmk(l1.size(),false);
  typedef Graphref_lst M_lst;
  sorted_sum_of_M_prod lMl;
  //  inc(tmk);
  while(!is_full(tmk)) 
    {
      inc(tmk);
      M_lst ml;
      int cnt = 0;
      BOOST_FOREACH(Graphref g,l1)
        {
          if(tmk[cnt])ml.push_back(g);
          cnt++;
        }
      bool zero_term = false;
      cout<<"Extracting terms"<<endl;
  //--------------------------------------------------------------------------------     
 cout<<"extraction procedure started for ml = "<<ml.size()<<" terms"<<endl;
  if(ml.size() > 1)
    {
      ///searchig for zero terms
      for (Graphref_lst::iterator it = ml.begin();it != ml.end();++it)
        for (Graphref_lst::iterator it2 = ml.begin();it2 != ml.end();++it2)
          if((it != it2)&&(isStrongIntersected(*it,*it2))) 
            {
              cout<<"ZERO TERM"<<endl;
              zero_term = true;
            }
      if(!zero_term)
        {
          /// sorting the list
          cout<<"Not sorted";
          BOOST_FOREACH(Graph g,ml)
            {
              cout<<" "<<num_vertices(g)<<" ";
            }
          cout<<endl;
          ml.sort(isSubgraph);
          //      ml.reverse();
          cout<<"After sorting";
          BOOST_FOREACH(Graph g3,ml)
            {
              cout<<" "<<num_vertices(g3)<<" ";
            }
          cout<<endl;
          lMl.push_back(ml);
         }
    }
  else
    {
      cout<<"Working with single term"<<endl;
      lMl.push_back(ml);
    }
  //--------------------------------------------------------------------------------     

      //   inc(tmk);
    }

  cout<<" number of permutations(lMl size) "<<lMl.size()<<endl;
  return lMl;
}


/// true if vertex have external leg
bool is_global_ext_vert(Graph::vertex_descriptor vdis,Graph& Gtmp)
{
  if(out_degree(vdis,Gtmp))
    {
      if(Gtmp.is_root())return(Gtmp[vdis].external);
      else  return(Gtmp.root()[Gtmp.local_to_global(vdis)].external);
    }
  else return false;
}

std::list<symbol> unique_symbols(ex in_ex)
{
  std::list<symbol> m_lst;
  for (const_preorder_iterator it = in_ex.preorder_begin();it!=in_ex.preorder_end();++it)
    if(is_a<symbol>(*it))      
      m_lst.push_back(ex_to<symbol>(*it));
  m_lst.unique();
  return m_lst;
}

// change unique rl in expression
ex fix_rl(ex in_ex,unsigned char newrl)
{
  using cnst::D;
  using cnst::C;
  using cnst::ep;
 
  unsigned char rl = clifford_max_label(in_ex);
  // replace dirac_ONE
  exset repls_one;
  if( in_ex.find(dirac_ONE(),repls_one)) 
    {
      cout<<"repls one  for rl "<<int(rl)<<"  "<<repls_one<<endl;     
      in_ex = in_ex.subs(dirac_ONE(rl)==dirac_ONE(newrl));
    }
  
  // replace dirac_slash
  exset repls_slash;
  if( in_ex.find(dirac_slash(wild(1),4-2*ep,rl),repls_slash)) // !!! CORRECT 4 to D dimensions
    {
      cout<<"repls slash  for rl "<<int(rl)<<"  "<<repls_slash<<endl;     
      exmap rl_subsmap;
      for(exset::const_iterator it = repls_slash.begin();it!=repls_slash.end();++it)
        {
          exmap rmap;
          it->match(dirac_slash(wild(3),4-2*ep,rl),rmap);
          ex tmp_mom = wild(3).subs(rmap);
          cout<<"repls "<<rmap<<" subs "<<(dirac_slash(tmp_mom,4-2*ep,newrl))<<endl;//.subs(rmap))<<endl;

          rl_subsmap[*it] = (dirac_slash(tmp_mom,4-2*ep,newrl));//.subs(rmap));
        }
      cout<<"repls final subs "<<rl_subsmap<<endl;
      in_ex = in_ex.subs(rl_subsmap);

    }

  // replace contrvariant dirac_gamma
  exset repls_gamma_contr;
  if( in_ex.find(dirac_gamma(varidx(wild(3),4-2*ep),rl),repls_gamma_contr))
    {
      cout<<"repls contr gamma  for rl "<<int(rl)<<"  "<<repls_gamma_contr<<endl;     
      exmap rl_subsmap;
      for(exset::const_iterator it = repls_gamma_contr.begin();it!=repls_gamma_contr.end();++it)
        {
          exmap rmap;
          it->match(dirac_gamma(varidx(wild(3),4-2*ep),rl),rmap);
          ex tmp_mom = wild(3).subs(rmap);
          cout<<"repls "<<rmap<<" subs "<<(dirac_gamma(varidx(tmp_mom,4-2*ep),newrl))<<endl;//.subs(rmap))<<endl;

          rl_subsmap[*it] = (dirac_gamma(varidx(tmp_mom,4-2*ep),newrl));//.subs(rmap));
        }
      cout<<"repls final subs "<<rl_subsmap<<endl;
      in_ex = in_ex.subs(rl_subsmap);

    }

  // replace covariant dirac_gamma
  exset repls_gamma_cov;
  if( in_ex.find(dirac_gamma(varidx(wild(3),4-2*ep,true),rl),repls_gamma_cov))
    {
      cout<<"repls co gamma  for rl "<<int(rl)<<"  "<<repls_gamma_cov<<endl;     
      exmap rl_subsmap;
      for(exset::const_iterator it = repls_gamma_cov.begin();it!=repls_gamma_cov.end();++it)
        {
          exmap rmap;
          it->match(dirac_gamma(varidx(wild(3),4-2*ep,true),rl),rmap);
          ex tmp_mom = wild(3).subs(rmap);
          cout<<"repls "<<rmap<<" subs "<<(dirac_gamma(varidx(tmp_mom,4-2*ep,true),newrl))<<endl;//.subs(rmap))<<endl;

          rl_subsmap[*it] = (dirac_gamma(varidx(tmp_mom,4-2*ep,true),newrl));//.subs(rmap));
        }
      cout<<"repls final subs "<<rl_subsmap<<endl;
      in_ex = in_ex.subs(rl_subsmap);

    }
      
  // replace dirac_gamma5
  exset repls_gamma5;
  if( in_ex.find(dirac_gamma5(rl),repls_gamma5))
    {
      cout<<"repls gamma5  for rl "<<int(rl)<<"  "<<repls_gamma5<<endl;     
    }

  // replace dirac_gammaR
  exset repls_gammaR;
  if( in_ex.find(dirac_gammaR(rl),repls_gammaR))
    {
      cout<<"repls gammaR  for rl "<<int(rl)<<"  "<<repls_gammaR<<endl;     
    }

  // replace dirac_gammaL
  exset repls_gammaL;
  if( in_ex.find(dirac_gammaL(rl),repls_gammaL))
    {
      cout<<"repls gammaL  for rl "<<int(rl)<<"  "<<repls_gammaL<<endl;     
    }
  return in_ex;
}

// fix representation labels(ferion line indices) in graph expressions
Graph fix_rl(Graph G_in)
{
  Graph G(G_in);
  std::list<rllist> rllst = get_rl_lst(G);
  Graph::edge_iterator ei,ei_end;
  Graph::vertex_iterator vi,vi_end;
  for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
    {
      rlchar local_rl = clifford_max_label(G[*ei].gexpr);
      int rl_offset = 0;
      int the_offset = 0;
      cout<<" where to find "<<rllst<<endl;
      BOOST_FOREACH(rllist l5,rllst)
        {
          rl_offset++;
          if(find(l5.begin(),l5.end(),local_rl)!=l5.end()) the_offset = rl_offset;
        }
      //      BOOST_ASSERT_MSG(the_offset>0," rl not found in propagator expression");
      if(the_offset >0)
        {
      G[*ei].gexpr = fix_rl(G[*ei].gexpr,40+the_offset);
      cout<<"fix "<< G[*ei].gexpr<<"    "<<fix_rl(G[*ei].gexpr,67)<<" rl: "<<int(clifford_max_label(G[*ei].gexpr))<<endl;
        }
    }
  for(tie(vi,vi_end) = vertices(G);vi!=vi_end;++vi)
    {
      rlchar local_rl = clifford_max_label(G[*vi].gexpr);
      int rl_offset = 0;
      int the_offset = 0;
      BOOST_FOREACH(rllist l5,rllst)
        {
          rl_offset++;
          if(find(l5.begin(),l5.end(),local_rl)!=l5.end()) the_offset = rl_offset;
        }
      //      BOOST_ASSERT_MSG(the_offset>0," rl not found in vertex expression");
      if(the_offset >0)
        {   
          G[*vi].gexpr = fix_rl(G[*vi].gexpr,40+the_offset);
          cout<<"fix "<< G[*vi].gexpr<<"    "<<fix_rl(G[*vi].gexpr,67)<<" rl: "<<int(clifford_max_label(G[*vi].gexpr))<<endl;
        }

    }

  for(tie(vi,vi_end) = vertices(G);vi!=vi_end;++vi)
    cout<<"ver "<<G[*vi].gexpr<<endl;
  cout<<"ex test "<<ex_from_graph(G)<<endl;
  return G;
}
// find set of representation labels with same spinor number
std::list<rllist> get_rl_lst(Graph G)
{
///////////////////////////////////////////////
  ///  _____________________________________  ///
  /// |                                     | ///
  /// | Expression from graph, not subgraph | ///
  /// |_____________________________________| ///
  ///                                         ///
  ///////////////////////////////////////////////

  if(G.is_root())
    {
      std::list<rllist> out_lst;
      GiNaC::numeric one1 = 1;
      ex out_expr = one1;
      cout<<"cliff "<<int(clifford_max_label(ex_from_graph(G)))<<endl;
      // add out spinor lines and vertices with external legs 
      Graph::vertex_iterator vi,vi_end;

      // mark vertices and edges as visited
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) G[*vi].visited = false;
      Graph::edge_iterator ei,ei_end;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei)
        {
          G[*ei].visited = false;
          //    cout<<"to "<<G[*ei].to_vertex<<" target"<<target(*ei,G)<<endl;
        }

      Graph::vertex_descriptor vdis;
      Graph::edge_descriptor edis;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && is_ext_spin_end(vertex(G[*ei].to_vertex,G),G) && !G[*ei].visited)
          {
            cout<<"rl found"<<endl;
            // list of rl_s for one spinor line
            rllist rl_lst;
            rl_lst.push_back(clifford_max_label(G[vertex(G[*ei].to_vertex,G)].gexpr));
            //    cout<<"01Vert in linesTesting out_expr:"<< out_expr<<endl;
            G[vertex(G[*ei].to_vertex,G)].visited = true;
            rl_lst.push_back(clifford_max_label(G[*ei].gexpr));
            // cout<<"01Edges  in lines Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
            edis = *ei;
            bool visited_spinor_line = false;
            while(!visited_spinor_line)
              {
                rl_lst.push_back(clifford_max_label(G[vertex(G[edis].from_vertex,G)].gexpr));
                //  cout<<"Vert in lines Testing out_expr:"<< out_expr<<endl;
                G[vertex(G[edis].from_vertex,G)].visited = true;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                bool lng = false;
                for(tie(out_ei,out_ei_end) = out_edges(vertex(G[edis].from_vertex,G),G);out_ei != out_ei_end;++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited)) 
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }

                if(lng == true)
                  {
                    rl_lst.push_back(clifford_max_label(G[*out_ei_tmp].gexpr));
                    //  cout<<"Edges in lines Testing out_expr:"<< out_expr<<endl;
                    G[*out_ei_tmp].visited = true;
                    edis = *out_ei_tmp;
                  }
                else visited_spinor_line = true;
              }
            rl_lst.unique();
            out_lst.push_back(rl_lst);
          }
      // spinor edges in loops
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && !G[*ei].visited)
          {
            rllist rl_lst;
            bool visited_loop = false;
            edis = *ei;
            while(!visited_loop)
              {
                rl_lst.push_back(clifford_max_label( G[vertex(G[edis].to_vertex,G)].gexpr));
                //cout<<"Vertieces in loops Testing out_expr:"<< out_expr<<endl;
                G[vertex(G[edis].to_vertex,G)].visited = true;
                rl_lst.push_back(clifford_max_label(G[edis].gexpr));
                // cout<<"Edges inloops Testing out_expr:"<< out_expr<<endl;
                G[edis].visited = true;

                bool lng = false;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                for(tie(out_ei,out_ei_end) = out_edges(vertex(G[edis].from_vertex,G),G); out_ei != out_ei_end; ++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited))
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)edis = *out_ei_tmp;		
                else visited_loop = true;
	      
              }
            rl_lst.unique();
            out_lst.push_back(rl_lst);
          }
      
      //cout<<"rl extracted from GRAPH "<<out_lst<<endl;
      //      cout<<"Not simplified: "<<out_expr<<endl;
      //cout<<"N = "<<out_expr.numer()<<endl;
      //cout<<"D = "<<out_expr.denom()<<endl;
      return out_lst;
    }
  ///////////////////////////////////////////////
  ///  _____________________________________  ///
  /// |                                     | ///
  /// | Expression from SUBgraph, not graph | ///
  /// |_____________________________________| ///
  ///                                         ///
  ///////////////////////////////////////////////
  else
    {
      //      print(G);
      GiNaC::numeric one1 = 1;
      ex out_expr = one1;
      cout<<"subgr Testing out_expr:"<< out_expr<<endl;
      // add out spinor lines and vertices with external legs 
      Graph::vertex_iterator vi,vi_end;
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) G[*vi].visited = false;
      Graph::edge_iterator ei,ei_end;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) G[*ei].visited = false;
        print(G);
        cout<<"After printing"<<endl;
      Graph::vertex_descriptor vdis;
      Graph::edge_descriptor edis;
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && is_ext_spin_end(target(*ei,G),G) && !G[*ei].visited)
          {
            out_expr *= G[target(*ei,G)].gexpr;
             cout<<"01Vert in linesTesting out_expr:"<< out_expr<<endl;
             G[target(*ei,G)].visited = true;
            out_expr *= G[*ei].gexpr;
             cout<<"01Edges  in lines Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
            edis = *ei;
            bool visited_spinor_line = false;
            while(!visited_spinor_line)
              {
                out_expr *= G[source(edis,G)].gexpr;
                  cout<<"Vert in lines Testing out_expr:"<< out_expr<<endl;
                G[source(edis,G)].visited = true;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                bool lng = false;
                for(tie(out_ei,out_ei_end) = out_edges(source(edis,G),G);out_ei != out_ei_end;++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited)) 
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)
                  {
                    out_expr *= G[*out_ei_tmp].gexpr;
                      cout<<"Edges in lines Testing out_expr:"<< out_expr<<endl;
                    G[*out_ei_tmp].visited = true;
                    edis = *out_ei_tmp;
                  }
                else visited_spinor_line = true;
              }
          }
      // spinor edges in loops
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if((G[*ei].type == 1) && !G[*ei].visited)
          {
            bool visited_loop = false;
            edis = *ei;
            while(!visited_loop)
              {
                out_expr *= G[target(edis,G)].gexpr;
                //   cout<<"Vertieces in loops Testing out_expr:"<< out_expr<<endl;
                G[target(edis,G)].visited = true;
                out_expr *= G[edis].gexpr;
                //   cout<<"Edges inloops Testing out_expr:"<< out_expr<<endl;
                G[edis].visited = true;

                bool lng = false;
                Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
                for(tie(out_ei,out_ei_end) = out_edges(source(edis,G),G); out_ei != out_ei_end; ++out_ei)
                  if((G[*out_ei].type == 1) && (!G[*out_ei].visited))
                    {
                      out_ei_tmp = out_ei;
                      lng = true;
                    }
                if(lng == true)edis = *out_ei_tmp;		
                else visited_loop = true;
	      
              }
          }
      // other edges and other verticies
      for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
        if(!G[*ei].visited)
          {
            out_expr *= G[*ei].gexpr;
            //   cout<<"Other_edges Testing out_expr:"<< out_expr<<endl;
            G[*ei].visited = true;
          }
      for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        if(!G[*vi].visited)
          {
            out_expr *= G[*vi].gexpr;
            cout<<*vi<< G[*vi].gexpr<< " Other Vertices Testing out_expr:"<< out_expr<<endl;
            G[*vi].visited = true;
          }

      cout<<"Expr extracted from GRAPH no."<<get_property(G,graph_num)<<endl;
      //      cout<<"Not simplified: "<<out_expr<<endl;
      //cout<<"N = "<<out_expr.numer()<<endl;
      //cout<<"D = "<<out_expr.denom()<<endl;
      //return lst(out_expr);

    }  
}

 std::ostream & operator<<(std::ostream & os, const expair& e)
 {
     os<<"{"<<e.rest<<","<<e.coeff<<"}";
     return os;
 }



// template <> 

std::ostream & operator<<(std::ostream & os, const rllist& e)
 {
   // print_context *p = get_print_context(os);
     rllist::const_iterator i = e.begin();
     rllist::const_iterator vend = e.end();
 
     if (i==vend) {
         os << "()";
         return os;
     }
 
     os << "(";
     while (true) {
       os<<int(*(i));
         ++i;
         if (i==vend)
             break;
         os << ",";
     }
     os << ")";
 
     return os;
 }


void print_rl_tree(ex in_ex)
{
  cout<<"cl_tr begin"<<endl;
  for(const_preorder_iterator it = in_ex.preorder_begin();it!= in_ex.preorder_end();++it)
    {
      if(is_a<clifford>(*it))
        {
          if(is_a<diracone>(*it))cout<<"dione!!! ";
          cout<<(*it)<<"cl_tr: "<<ex_to<clifford>(*it).get_representation_label()<<endl;
        }
    }
}

// from subs map like 
//    k1+l1=p1
//    k2+l1=p3
//    get map for subs only for external momentums
exmap map_to_subs(exmap in_map)
{
  //symbol l1("l1"),l2("l2"),l3("l3"),l4("l4"),l5("l5");
  //const ex lm[] = {l1,l2,l3,l4,l5};
  int n = 0;
  //std::list<ex> m_lst_lst;
  //list of all momentum symbols in external lines 
  std::set<symbol,ex_is_less> m_set;
  //   BOOST_FOREACH(ex e478,ext_vec)
  for(exmap::iterator map_it = in_map.begin();map_it!=in_map.end();++map_it)  
   {
     ex e478 = map_it->first;
     cout<<"mom of vert: "<<endl;
     std::copy(e478.preorder_begin(),e478.preorder_end(),std::ostream_iterator<ex>(cout, "\n"));
     cout<<"my way:"<<endl;
      
      for (const_preorder_iterator it = e478.preorder_begin();it!=e478.preorder_end();++it)
        {
          if(is_a<symbol>(*it))
            if(cnst::km.has(ex_to<symbol>(*it))) // insert only loop momentums
              {
                m_set.insert(ex_to<symbol>(*it));
                cout<<"---> "<<(*it)<<" is a symbol: "<<is_a<symbol>(*it)<<endl;
              }
        }
        
          //      cout<<"solution subs "<<(*it)<<" = "<<lsolve(e478==lm[n],(*it))<<endl;
    }
  cout<<"not unique symbols: "<<m_set.size()<<endl;
  /*  BOOST_FOREACH(symbol e634,m_set)
    {
      cout<<"not unique "<<e634<<endl;
      }*/
  //  m_set.sort();
  // m_set.unique();
  BOOST_FOREACH(symbol e634,m_set)
    {
      cout<<"unique "<<e634<<endl;
    }
  /// creation of momentum index map
  //std::map<symbol,int> mom_idx;
  exhashmap<int> mom_idx;
  BOOST_FOREACH(ex e789,m_set)
    {
      //      BOOST_FOREACH(ex e45,ext_vec)
      for(exmap::iterator map_it = in_map.begin();map_it!=in_map.end();++map_it)  
        {
          ex e45 = map_it->first;
          
          if(e45.has(e789))
            {
              cout<<e45<<" has "<<e789<<endl;
              mom_idx[ex_to<symbol>(e789)]++;
            }
        }
    }
  for(exhashmap<int>::iterator it = mom_idx.begin();it != mom_idx.end();++it)
    cout<<"mom: "<<it->first<<" mom_idx: "<<it->second<<endl;
  /// finding of minimal momentum index in each ext_mom
  /// and mom map creation
  //list in form (ext_mom,min_symb_mom,idx),where min is a momentum with minimal index
  std::vector< tuple<ex,symbol,int> > em_m_lst;
  // BOOST_FOREACH(ex e45,ext_vec)
  
  for(exmap::iterator map_it = in_map.begin();map_it!=in_map.end();++map_it)  
    {
      ex e45 = map_it->first;
  
      symbol buf_mom;
      int buf_idx = 1000;
      BOOST_FOREACH(symbol e789,m_set)
        {
          if(e45.has(e789) && (mom_idx[e789] < buf_idx))
            {
              buf_idx = mom_idx[e789];
              buf_mom = e789;
            }
        }
      em_m_lst.push_back(make_tuple(e45,buf_mom,buf_idx));
    }
  sort(em_m_lst.begin(),em_m_lst.end(),nops_less_comp);
  sort(em_m_lst.begin(),em_m_lst.end(),idx_less_comp);

  typedef tuple<ex,symbol,int> mom_tuple;
  
  exmap mom_subs_map,mom_subs_back_map;    
  exvector l_mom_vec;
  enum{ext_mom_expression,mom_symbol,mom_index};
  ex buf_momentum = 0;
  BOOST_FOREACH(mom_tuple m_pair,em_m_lst)
    {
      //      if(n<(em_m_lst.size()-1))
        {
          mom_subs_map[m_pair.get<mom_symbol>()] = lsolve(m_pair.get<0>()==in_map[m_pair.get<0>()]/*lm[n]*/,m_pair.get<mom_symbol>());
          mom_subs_back_map[in_map[m_pair.get<0>()]] = m_pair.get<ext_mom_expression>();
          cout<<"Ext mom: "<<m_pair.get<0>()<<" min "<<m_pair.get<1>()<< " index "<<mom_idx[m_pair.get<1>()]<<endl;
          // l_mom_vec.push_back(lm[n]);
          //buf_momentum+=lm[n];
        }
        /*else
        {
          mom_subs_map[m_pair.get<mom_symbol>()] = lsolve(m_pair.get<0>()==-buf_momentum,m_pair.get<mom_symbol>());
          //          mom_subs_back_map[lm[n]] = m_pair.get<ext_mom_expression>();
          cout<<"Last mom: "<<m_pair.get<0>()<<" min "<<m_pair.get<1>()<< " index "<<mom_idx[m_pair.get<1>()]<<endl;
          l_mom_vec.push_back(-buf_momentum);
          }*/
      n++;
    }
  //  cout<<colors::blue<<"SMOM:  "<<in_ex.subs(mom_subs_map)<<colors::def_col<<endl;

  //cout<< "map 1 "<<P_map<< " map 2 "<<P_map.subs(P_map)<<endl;
  exmap tmp_map = mom_subs_map;
  for(exmap::iterator map_it = mom_subs_map.begin();map_it!=mom_subs_map.end();++map_it)  
    map_it->second = (map_it->second).subs(tmp_map);
  return mom_subs_map;//,mom_subs_back_map);
}

bool idx_less_comp(tuple<ex,symbol,int> left,tuple<ex,symbol,int> right)
{
  return(get<2>(left) < get<2>(right));
}

bool nops_less_comp(tuple<ex,symbol,int> left,tuple<ex,symbol,int> right)
{
  return(get<0>(left).nops() > get<0>(right).nops());
}



// return trace over all rl
ex trace_ng(ex in_ex)
{
  std::set<rlchar> rl_set;
 for(const_preorder_iterator it = in_ex.preorder_begin();it!= in_ex.preorder_end();++it)
    {
      if(is_a<clifford>(*it) )
        {
          rl_set.insert( ex_to<clifford>(*it).get_representation_label());
        }
    }
 BOOST_FOREACH(rlchar rl,rl_set)
   {
     cout<<"trace before "<<in_ex<<endl;
     in_ex = dirac_trace(in_ex,rl);
   }
 cout<<"trace final "<<in_ex<<endl;
   return in_ex;
}
// true if in_ex with rl has not dirac_ONE clifford elements 
bool has_not_dirac_ONE(ex in_ex,unsigned char rl)
{
  BOOST_ASSERT_MSG(!is_a<add>(in_ex),"Finding clifford units in ADD, not in MUL");
  int cntr = 0;
   for(const_preorder_iterator it = in_ex.preorder_begin();it!= in_ex.preorder_end();++it)
    {
      if(is_a<clifford>(*it) && (!is_a<diracone>(*it)))
        {
          if( (ex_to<clifford>(*it).get_representation_label()) == rl)          cntr++;
        }
    }
   return (cntr>0);
}

//remove dirac_ONE with rl, if exists non dirac_ONE terms with this rl
ex remove_ONE_in_clifford(ex in_ex,unsigned char rl)
{
  in_ex = in_ex.expand();
  ex tmp_ex = in_ex;
  if(is_a<add>(in_ex))
    {
      cout<<"ONE add " <<in_ex<<endl;
      for(const_iterator it = in_ex.begin();it!=in_ex.end();++it)
        {
          tmp_ex = tmp_ex.subs((*it) == remove_ONE_in_clifford(*it,rl));
        }
      return tmp_ex;
    }
  else if(is_a<mul>(in_ex))
    {
      cout<<"ONE mul  "<<has_not_dirac_ONE(in_ex,rl)<<"   " <<in_ex<<endl;
 exset repls;
      in_ex.find(dirac_ONE(rl),repls);
      cout<<"rl "<<rl<<" repls "<<repls<<endl;
      if(has_not_dirac_ONE(in_ex,rl)) return  in_ex.subs(dirac_ONE(rl) == 1);
     
    }
  else       cout<<"ONE x3 " <<in_ex<<endl;
  return in_ex;
}

ex trace_rl(ex in_ex,rllist l_in)
{
  try{
  using cnst::D;
  using cnst::C;
  using cnst::ep;
  const unsigned char newrl = 41;
  cout<< "repls in_ex: "<<in_ex<<endl<<endl;
  BOOST_FOREACH(rlchar rl,l_in)
    {

      // replace dirac_ONE
      exset repls_one;
      if( in_ex.find(dirac_ONE(rl),repls_one)) 
        {
          cout<<"repls one  for rl "<<int(rl)<<"  "<<repls_one<<endl;     
        }

      // replace dirac_slash
      exset repls_slash;
      if( in_ex.find(dirac_slash(wild(1),4-2*ep,rl),repls_slash)) // !!! CORRECT 4 to D dimensions
        {
          cout<<"repls slash  for rl "<<int(rl)<<"  "<<repls_slash<<endl;     
          exmap rl_subsmap;
          for(exset::const_iterator it = repls_slash.begin();it!=repls_slash.end();++it)
            {
              exmap rmap;
              it->match(dirac_slash(wild(3),4-2*ep,rl),rmap);
              ex tmp_mom = wild(3).subs(rmap);
              cout<<"repls "<<rmap<<" subs "<<(dirac_slash(tmp_mom,4-2*ep,newrl))<<endl;//.subs(rmap))<<endl;

              rl_subsmap[*it] = (dirac_slash(tmp_mom,4-2*ep,newrl));//.subs(rmap));
            }
          cout<<"repls final subs "<<rl_subsmap<<endl;
          in_ex = in_ex.subs(rl_subsmap);

        }

      // replace contrvariant dirac_gamma
      exset repls_gamma_contr;
      if( in_ex.find(dirac_gamma(varidx(wild(3),4-2*ep),rl),repls_gamma_contr))
        {
          cout<<"repls contr gamma  for rl "<<int(rl)<<"  "<<repls_gamma_contr<<endl;     
          exmap rl_subsmap;
          for(exset::const_iterator it = repls_gamma_contr.begin();it!=repls_gamma_contr.end();++it)
            {
              exmap rmap;
              it->match(dirac_gamma(varidx(wild(3),4-2*ep),rl),rmap);
              ex tmp_mom = wild(3).subs(rmap);
              cout<<"repls "<<rmap<<" subs "<<(dirac_gamma(varidx(tmp_mom,4-2*ep),newrl))<<endl;//.subs(rmap))<<endl;

              rl_subsmap[*it] = (dirac_gamma(varidx(tmp_mom,4-2*ep),newrl));//.subs(rmap));
            }
          cout<<"repls final subs "<<rl_subsmap<<endl;
          in_ex = in_ex.subs(rl_subsmap);

        }

      // replace covariant dirac_gamma
      exset repls_gamma_cov;
      if( in_ex.find(dirac_gamma(varidx(wild(3),4-2*ep,true),rl),repls_gamma_cov))
        {
          cout<<"repls co gamma  for rl "<<int(rl)<<"  "<<repls_gamma_cov<<endl;     
          exmap rl_subsmap;
          for(exset::const_iterator it = repls_gamma_cov.begin();it!=repls_gamma_cov.end();++it)
            {
              exmap rmap;
              it->match(dirac_gamma(varidx(wild(3),4-2*ep,true),rl),rmap);
              ex tmp_mom = wild(3).subs(rmap);
              cout<<"repls "<<rmap<<" subs "<<(dirac_gamma(varidx(tmp_mom,4-2*ep,true),newrl))<<endl;//.subs(rmap))<<endl;

              rl_subsmap[*it] = (dirac_gamma(varidx(tmp_mom,4-2*ep,true),newrl));//.subs(rmap));
            }
          cout<<"repls final subs "<<rl_subsmap<<endl;
          in_ex = in_ex.subs(rl_subsmap);

        }
      
      // replace dirac_gamma5
      exset repls_gamma5;
      if( in_ex.find(dirac_gamma5(rl),repls_gamma5))
        {
          cout<<"repls gamma5  for rl "<<int(rl)<<"  "<<repls_gamma5<<endl;     
        }

      // replace dirac_gammaR
      exset repls_gammaR;
      if( in_ex.find(dirac_gammaR(rl),repls_gammaR))
        {
          cout<<"repls gammaR  for rl "<<int(rl)<<"  "<<repls_gammaR<<endl;     
        }

      // replace dirac_gammaL
      exset repls_gammaL;
      if( in_ex.find(dirac_gammaL(rl),repls_gammaL))
        {
          cout<<"repls gammaL  for rl "<<int(rl)<<"  "<<repls_gammaL<<endl;     
        }


    }
  in_ex =   in_ex.simplify_indexed();
  
  //in_ex = remove_dirac_ONE(in_ex,0,1);
  //cout<<"repls wo trace "<<in_ex<<" "<<int(clifford_max_label(in_ex,true))<<endl;
  print_rl_tree(in_ex);
  //  in_ex = remove_dirac_ONE(in_ex,41,1);
  in_ex = canonicalize_clifford(in_ex);
  //in_ex = remove_ONE_in_clifford(in_ex,newrl);
  cout<<"repls exper "<<in_ex<<endl;
  in_ex = dirac_trace(in_ex,41);
  //in_ex =   in_ex.simplify_indexed();
  cout<<"repls max lab "<<in_ex<<"  "<<int(clifford_max_label(in_ex,true))<<endl;

  return in_ex;
 }
  catch(std::exception &p)
    {
      //     cout<<"******************************************************"<<endl;
      //cout<<">>>  ERROR in class AFDC:  "<<setw(40)<<std::internal<<p.what()<<endl;
      //cout<<"******************************************************"<<endl;
     throw std::logic_error(std::string("In function \"TRACE_RL\":\n |___> ")+p.what());
     }
}
